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Abstract 

We develop a method for approximating the Grobner basis of the ideal 
of polynomials which vanish at a finite set of points, when the coordinates 
of the points are known with only limited precision. The method con- 
sists of a preprocessing phase of the input points to mitigate the effects 
of the input data uncertainty, and of a new "numerical" version of the 
Buchberger-Moller algorithm to compute an approximation GB to the 
exact Grobner basis. This second part is based on a threshold-dependent 
procedure for analyzing from a numerical point of view the membership 
of a perturbed vector to a perturbed subspace. With a suitable choice 
of the threshold, the set GB turns out to be a good approximation to a 
"possible" exact Grobner basis or to a basis which is an "attractor" of 
the exact one. In addition, the polynomials of GB are "sufficiently near" 
to the polynomials of the extended basis, introduced by Stetter, but they 
present the advantage that LT(GB) coincides with the leading terms of a 
"possible" exact case. The set of the preprocessed points, approximation 
to the unknown exact points, is a pseudozero set for the polynomials of 
GB. 

Keywords: Ideal of points, Grobner basis, perturbed data, numerical algo- 
rithm. 



1 Introduction 

Let P = K[x\, . . . , x s ] be the polynomial ring in s indeterminates over a field 
K . In this paper we develop a method for approximating the Grobner basis of 
the ideal IcFof polynomials which vanish at a finite set of points, when the 
coordinates of the points are known with only limited precision. In particular 
we analyze the case K = M. 

In the exact case, that is when we deal with a set of unperturbed points V ', 
the problem has been deeply analyzed by several authors (see for example 2 , 
[S], [IH]) [H]) and the Grobner basis of the ideal T(V) of the points V can be 
computed, for example, by the Buchberger-Moller algorithm (pQ, [4], [10]). It is 
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well known that, given a term ordering, each set of m points V corresponds to 
a closed set TV = {ii, . . . , t m } such that each polynomial g of the Grobner basis 
GB oiJ(V) has the form 



with leading term t and suitable coefficients Cj. The set AT, called normal set, 
is the basis of the quotient ring P/I{V) and the leading term t belongs to the 
corner set C [N] . 

Nevertheless, it is also well known that the problem of computing the Grobner 
basis of an ideal of points is "discontinuous" ([TO], [H]), i.e. small perturbations 
of the points can cause structural changes in the basis, as illustrated in Example 
1.1. For this reason exact methods applied blindly to perturbed "real- world" 
data can produce meaningless results. 

Example 1.1. Given the term ordering DegLex and the exact points p\ — 
(1,1), P2 = (3,2) andp3 = (5, 3), the Grobner basis of the ideal Z of polynomials 
in R[x, y] vanishing at these points is 



This is an obvious result since the points V are not aligned and so GB struc- 
turally differs from GB. Nevertheless, since the points V are not aligned because 
of data errors, from the computational point of view it is more interesting to 
approximate the Grobner basis GB than to compute GB, given as input the set 
V. 

In the general case, if the input points arise from real-world data, their co- 
ordinates are approximations to the unknown exact values, and only an error 
estimation is known. For this reason, each point p of V represents a "cloud" 
of points. More precisely, given an input point p known with uncertainty So, 
each point p which differs from p componentwise by less than sq can be cho- 
sen as an input point computationally equivalent to p, and so it will be called 
"admissible input point" . Analogously, given a set of perturbed input points 
V, an "admissible input set " is a set consisting of admissible input points 
computationally equivalent to V . Obviously, given an error estimation sg, there 
exist several admissible input sets. These remarks combined with the struc- 
tural discontinuity show that the computation of the Grobner basis of ideals 
of perturbed points is a very tricky problem. In fact, since the exact points 





with normal set N = {1, y, y 2 }. 



If data errors perturb the point p% so that the set of input points is V 
{Pii Vii (5-1,3)}, then the Grobner basis GB of the ideal T(V) is: 




with normal set Af — {1, y, x}. 
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are unknown and the Grobner basis can change choosing different admissible 
input sets, it can be difficult to identify the case to be approximated. For this 
reason, we define some criteria to choose the "reference" Grobner basis to ap- 
proximate, pointing out its structure, that is its leading terms and the supports 
of its polynomials. 

We define "reference basis" a Grobner basis GB generating an ideal of 
polynomials vanishing at some admissible input set, whose structure represents 
all the admissible input sets. Since the structure of a Grobner basis and the 
elements of its normal set are closely connected, we require that the normal 
set of a reference basis can be associated to all the admissible input sets, that 
is, analogously to the exact case, we require that its terms provide independent 
vectors if evaluated on the points of any admissible input set. A normal set with 
this property will be called "reference normal set" . Note that a reference 
basis could coincide with the exact basis, since the set of exact points belongs 
to the cloud of the admissible input sets. 

In order to approximate a reference basis, our method computes a set of 
polynomials GB, whose supports belong to A7"u CfA/], where AT is a reference 
normal set. In some particular case, GB is a reference Grobner basis, and, 
in general, it is a good approximation to a reference basis. In fact, we show 
that, in the general case, GB = {g l7 . . . ,g~ n } structurally corresponds to a 
reference basis GB = {gi, . . . ,g n }, that is there is an one-to-one correspondence 
between GB and GB which preserves the leading terms and the supports of 
each polynomial. In addition, for each i — l...n, 'g i G GB and gi G GB 
have "similar" coefficients and we estimate the difference between 'g i and gi. 
Moreover, we show that the elements of GB are "near" to the polynomials of 
the extended basis GBe introduced by Stetter in [12] , 

Unfortunately, some "degenerate" cases can occur, when there are no ref- 
erence bases, that is when all the Grobner bases, with reference normal set, 
are associated to points which are not admissible input points. Nevertheless, in 
these degenerate cases the numerical tests point out that also the normal set 
of the exact points cannot represent, from a numerical point of view, all the 
admissible input sets. In addition, the tests suggest that the computed normal 
set 77 seems to be the numerically stable normal set "nearest" to the exact 
one. The analyzed degenerate examples are not presented in this paper since 
the degenerate cases will be studied in a future work, whereas in this article we 
shall limit ourselves to "non-degenerate" cases. 

Our method consists of two parts. First of all, the input points are prepro- 
cessed, in order to mitigate some negative effects of the data errors. The possible 
"splittings" of the coordinates due to the errors are removed and the prepro- 
cessing algorithm computes a new admissible input set V , by replacing with a 
unique suitable value the perturbed coordinates which differ by less than the 
error estimation. The preprocessed points turn out to be a set of pseudozeros 
[T3] for the polynomials of GB. 

The second part of our algorithm constitutes the main innovative part of 
this paper. It is a numerical version of the Buchberger-Moller algorithm which 
computes a set of polynomials GB, approximation to a reference basis. The 
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numerical Buchberger-Moller algorithm is based on the following notion of "ap- 
proximate" linear dependence of vectors. From the numerical point of view, 
a vector v depends on the vectors {vi, . ..,«*} if it is "sufficiently" near to 
an element of the subspace Spanjui, . . . , Vk}- Since this "numerical" depen- 
dence of vectors is defined choosing a threshold e, the numerical version of the 
Buchberger-Moller algorithm is threshold-dependent. We have to choose a suit- 
able value of the threshold for computing a good approximation to a reference 
basis. The threshold choice is one of the trickier aspects of the implementation 
of our method. 

This paper is organized as follows. In Section 2 we describe the threshold- 
dependent approximation algorithm and in Section 3 we present an heuristic 
approach for the choice of this threshold, whereas a more detailed analysis is 
left to future works. The approximation properties of the computed set GB 
are analyzed in Section 4, and a comparison with some results presented by 
Stetter is described in Section 5. In Section 6 we report two examples in order 
to illustrate and clarify explicitly the behaviour of our method. Appendix A 
contains a description of the strategy for preprocessing the input points used 
in our examples and Appendix B presents two error upper bounds useful in the 
estimation of the goodness of the basis approximation. 

Notation. Later on "G-basis" means "Grobner basis" . We denote the points 
with the letter p and the tuple of points with V = {pi,P2 ■ ■ •}■ If g{x\, ■ ■ ■ , x s ) 
is a polynomial in R[a;i, . . . , x s ], g(p) is the evaluation of g at p and g(V) is the 
vector whose components are the values of g on the elements of V . In addition 
we denote with ||^||2 the 2-norm of a vector v. 

2 The approximation method 

In this section we describe our method for approximating a reference basis with 
a set of polynomials GB, when we deal with perturbed data in non-degenerate 
cases, that is when a reference basis exists. 

Our method computes the set of polynomials GB, working on a set of points 
obtained after preprocessing the perturbed input data. The goal of the prepro- 
cessing phase is to minimize the damage caused by possible "splittings" of the 
data due to the input perturbations, replacing by a unique value the coordi- 
nates which differ by less than the error estimation, since they are coincident 
from the computational point of view. The preprocessed set V depends on the 
preprocessing technique; the strategy used in our tests is described in Appendix 
A. 

Even if the exact Buchberger-Moller (BM) algorithm 4 applied to the pre- 
processed data sometimes computes a good approximation to a reference basis, 
this is not true in general. The simple preprocessing of the input data is not 
sufficient to obtain a numerically significant set of polynomials, and so it is 
necessary to develop a new approximation algorithm. 
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It is well known that, given a term ordering a and a set of m points 
V = {pi, . . . ,p m }, not perturbed by errors, the BM algorithm computes the 
er-Grobner basis GB of the ideal of the polynomials which vanish at the set V . 
At each step, after choosing the term t, the BM algorithm checks if the vector 
t(V) is linearly independent of the vectors {ti(V), . . . , ifc('P)}, where {t%, . . . , £&} 
is the subset of the normal set computed at the previous steps, and t > a ti, 
i = 1, . . . , k. If this is the case, t is a new term of the normal set. Other- 
wise, the BM algorithm builds a new polynomial g of the G-basis GB such that 
g = t — J2i=i c rfi- The coefficients a, i — 1, . . . , k, are the components of the 
vector c that satisfies 

M k {V)c = t{V), (1) 

where M^iV") is the matrix [ti(V), ■ . ■ , tk(V)], whose columns are linearly inde- 
pendent, by construction. 

The above test of linear dependence is crucially affected by even very small 
variations of the input data. Therefore, when we deal with approximate input 
points, different choices of admissible data may lead to different choices in the 
BM algorithm. Nevertheless, if the data are affected by errors, it is possible 
to check the linear dependence of perturbed vectors from a numerical point of 
view. Intuitively, a perturbed vector v belongs, in some approximate sense, to 
a subspace W, spanned by perturbed vectors, if v is "sufficiently near" to an 
element of W. This idea of "numerical" membership of a vector to a subspace, 
formalized and analyzed by the author in [TJ, is based on the following definition. 



Definition 2.1 Given a subspace W and a threshold e > 0, a vector v "numer- 
ically" lies in W if 

tk <e 

N| 2 " ' 

where r is the component of v orthogonal to W. In this case we write v S € W. 

By exploiting this concept of "numerical" membership of a subspace, we 
develop the numerical Approximate Buchberger Moller algorithm (denoted by 
ABM) for computing the approximation GB to a reference basis. Working on 
the preprocessed input points P, the ABM algorithm checks, at each step, if a 
perturbed vector lies numerically in a subspace spanned by a perturbed basis. 
If the vector "numerically" lies in the subspace, then the algorithm builds a new 
polynomial of GB, even if the vector does not exactly belong to the subspace. 
Otherwise the ABM algorithm inserts a new element into the normal set 77. 

The ABM algorithm can be described as follows. We check the numerical 
membership of a vector v to a subspace W — Span{iui, . . . , Wk} computing the 
residual r S W 1 - of the least square problem [wi ... Wk]x = v, since, in this way, 
we can more easily estimate the difference between GB and a reference basis, 
as shown in Section 4. 
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The ABM Algorithm. 

Input. The term ordering cr, the preprocessed input points V — {p 11 . . . ,p m } 
and the threshold e. 

Output. The set of polynomials GB and the set of power products J7. 

Let t\ , . . . , tk be the terms of 77 computed at the previous steps and let t be the 
current power product, t > a t\, . ■ ■ , tk, chosen analogously to the BM algorithm. 

1. Given the matrix Mk{V) — pi (V), ... ,tk{V)], solve the least square prob- 
lem 

M k (V)c = t(V). (2) 

2. Let r = t(V) — M^c be the residual of the problem @, then 

• if ||f||2/PCP)||2 > e , put the term t into the set 77; 

• otherwise, put the polynomial ~g = t — Y]- Citi into GB. <^> 

Note that, analogously to the exact case, the supports of the polynomials of 
GB consist of terms of AfLiC[N~\ and the leading terms belong to C[Sf\. 

Obviously, the output sets GB and M of the ABM algorithm depend on the 
value of e and we have to choose a suitable value of the threshold so that GB 
is a good approximation to a reference basis. A detailed analysis of the choice 
of suitable thresholds will be presented in a future work. In the next section we 
illustrate an heuristic approach based on the analysis of some properties of M 
and GB. 

3 The threshold choice 

Intuitively, if we choose too small a value of e, the ABM algorithm might not 
recognize numerically dependent vectors and so the set GB could be very sensi- 
tive to data perturbations. Vice versa, if the value of the threshold is too large 
a vector may be considered as numerically belonging to a subspace even if it is 
too "far" from it. In this section we present an heuristic strategy for choosing a 
suitable value of the threshold analyzing the computed output sets. If the value 
of e does not produce satisfactory output sets, we repeat the procedure with a 
new threshold. 

As mentioned earlier, GB is a good approximation to a reference basis GB 
if each polynomial of GB corresponds to a polynomial of GB with the same 
support and similar coefficients and vice versa. In order to obtain the same 
structure of a reference basis, since we have excluded degenerate cases, we re- 
quire that Af is the normal set of a reference basis, and so we have to choose a 
value of e such that the following properties hold. 
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PI. Wc require that 77 has m elements. In this case, since the ABM algorithm 
analyzes the power products using the same strategy of the exact BM 
algorithm, 77 = {ti, . . . ,t m } turns out to be, by construction, a closed 
set. For this reason, as is well known [17], Af is the basis of the quotient 
ring P/I(fP) for each polynomial ideal T(fP) with m simple zeros V such 
that {ii(V), . . ■t m (V)} are linearly independent vectors, that is N is a 
normal set associated to such set V. 

P2. We require that for each term t the relative residual ||r||2/PCP)||2 of the 
problem with right side t(V) must be sufficiently greater or less than 
the threshold, that is well separated from it. This condition guarantees 
that small data perturbations do not affect each residual sufficiently to 
change its position with respect to e and so 77 turns out to be invariant 
for admissible input sets. It follows that the elements of Af provide inde- 
pendent vectors if evaluated on admissible input sets, and so, if PI holds, 
77 is a normal set for the admissible input sets, i.e. it is a reference normal 
set. 

P3. We require that, given the term ordering a, 77 is the normal set of a 
cr-Grobner basis GB of an ideal of admissible input points, that is 77 is 
the normal set of the reference G-basis GB. 

The properties PI and P2 can be directly checked on the set Af and on the 
computed relative residuals, whereas it can be more difficult to investigate if AT 
satisfies P3, since several cases can occur. 

If the following system, consisting of the polynomials of GB, 

9i = 

: (3) 

9n = 0, 

vanishes at V then GB is the cr-Grobner basis of the ideal I(V) with reference 
normal set Af, since in this case the ABM and BM algorithms have the same 
behaviour. Since V is an admissible input set, then GB is a reference basis and 
the property P3 is obviously satisfied. 

Otherwise we can check the property P3 using the following intuitive idea. 
The value of the threshold is suitably chosen if GB approximates the reference 
basis GB — {g\, . . . ,g n }, with normal set Af, corresponding to a set of points 
V near to V, that is if the function: 

^:p~ &&>)>••• »3„(P)), ( 4 ) 
is a good approximation to the function: 

F ip>-> (gi(p), ■ ■ -,9n(p)), 

which vanishes at V . For this reason, if the threshold is suitably chosen, we 
can suppose that Newton method [6] works analogously if applied to the poly- 
nomial system F(p) = or to F(p) = 0. Since Newton iteration applied to 
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the equation F(p) = computes descent directions for ||-F(p)||2 we can assume 
that it computes descent directions for ||F(p)||2 too. For this reason the set 
V + , computed by Newton iteration working on F with initial points V, can 
be considered a good approximation to V . It follows that the property P3 is 
"numerically" verified if the set V + is sufficiently "near" to V. 

Since the properties P1-P3 must be checked on the computed sets M and 
GB, we need a value of the threshold to begin the computation. The first value 
of the threshold can be chosen as follows. We apply the exact BM algorithm to 
the preprocessed points V and we compute for each term t of the exact normal 
set No the relative residual of the least square problem with right side t(P). We 
choose the first value of e in order to eliminate from J7o, by means the ABM 
algorithm, the power products with too small relative residuals. 

The following example illustrates the effects of different thresholds. 

Example 3.1 

Let V = {(1, 1), (3, 2), (5.1, 3)} be the perturbed points of Example 
1.1. Since the preprocessing phase does not change the perturbed points, we 
have V — V, corresponding to the exact normal set A/"o = {1, y, %}• Since the 
relative residuals of the terms y and x computed on V are equal to 0.37 and 
0.0068, respectively, we choose e > 0.0068, to obtain a set GB different from 
the exact G- basis of X('P). 

Choosing 0.0068 < e < 0.37, the polynomial x — 2.05y + 1.06 is inserted into 
GB and the term y 2 is analyzed. Since the relative residual of y 2 is equal to 
0.0824 we have two possibilities. 

If 0.0824 < e < 0.37, the nonzero relative residual of y 2 is treated as if it 
were equal to zero and the following set of polynomials is computed: 



GB 1 = 



x- 2.05y+ 1.06 
y 2 - Ay + 3.3. 



Since GB X vanishes at v\ = {(4.7072, 2.8165), (1.3595, 1.1835)} which differs 
in cardinality from V, we conclude that GB\ is not a good approximation to 
GB and that the threshold has not been suitably chosen. 

If 0.0068 < e < 0.0824, then the polynomial y 2 - Ay + 3.3 does not belong 
to GBi and the term y 2 is a new element of the normal set. So in this case the 
ABM algorithm computes the set GB2: 



GB 2 



a; -2.05y+ 1.06 
y 3 -6y 2 + lly-6 



G-basis of the ideal of points v\ = {(0.983, 1), (3.03, 2), (5.083, 3)}, very 
near to the set of points V, and so we consider the threshold suitably chosen. 
Note that GBi is very similar to the G-basis of the ideal of the exact points of 
Example 1.1. <> 

Remark. If the ABM algorithm cannot compute normal sets which satisfy 
properties P1-P3, then we are dealing with "degenerate" cases. The numerical 



8 



tests suggest that, in the degenerate examples, the relative least square residuals 
for the exact normal sets associated to all the admissible input points are small. 
Unfortunately, since it is rather difficult, if not impossible, to distinguish the case 
where the small residuals computed on the preprocessed points are caused by 
data perturbations from the case where they are related to degenerate problems, 
we cannot check, a priori, if we are dealing with degenerate cases. 



4 Comparison with a reference basis 

In this section we estimate the difference between a reference basis and the set 
GB computed, given the term ordering a, by the ABM algorithm working on 
the preprocessed points V . We assume that the value of the threshold e has been 
suitably chosen and that "non-degenerate" problems are analyzed, that is the 
computed set AT is a reference normal set and there exists a reference cr-Grdbncr 
basis GB with normal set N. We show that GB is a good approximation to a 
reference basis GB. 

First of all, note that if at each step of the ABM algorithm the analyzed 
relative residual is equal to zero or greater than the threshold, then the ABM 
algorithm works on the preprocessed points V like the BM algorithm. In this 
case the set GB is the exact G-basis of the ideal I(P) and so, since V is an 
admissible input set, GB coincides with a reference basis. 

In the general case, it can happen that, at some step of the ABM algorithm, 
the relative residual associated to t is greater than zero and less than or equal 
to e. In this case the polynomial g = t — y\. —1 Cjij is the new element of the set 
GB, even though it does not vanish at the points V , since 

k 

g{V) = t{V)-Y J cMT 5 )=r, 
i=l 

where f is the nonzero residual of ([2]), and so GB is not the G-basis of the 
ideal T(V). Nevertheless, since AT is also the normal set of a reference basis 
GB, GB structurally corresponds to GB. In addition, in the following we 
prove that the polynomials of GB and GB have similar coefficients. For these 
reasons, as mentioned earlier, the computed set GB can be considered a good 
approximation to a reference basis. 

Let V be the admissible input points which are zeros of the polynomials of 
a reference basis GB. By construction, if g S GB corresponds to g £ GB, the 
vectors c and c of the coefficients of g and g satisfy respectively 

M k {V)c = t(V)-r and M k (V)c = t(V), 

where r is the residual of the least square problem ([2]) using the preprocessed 
points. From the analysis of the sensitivity of a least square problem ([BJ, [5]), 
since the residual of the system ([1]) is zero, we have 

gg#> [„ + «,], (5 ) 



|c|| 2 ~ l-e M K 2 (M k {V))' 
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where 



e* = 



£m = 



\\M k {V)-M k (V)\\ 2 
\\Mk(V)\\ 2 



and K2(Mf.(V)) is the 2-condition number of the matrix MkiV) [3J. It is well 
known 3: that, in general, the upper bound ([5]) overestimates the sensitivity of 
the linear systems. Nevertheless, the formula ([5]) suggests that e t and ejj can be 
interpreted as a measure of goodness of approximation. If 5r is the maximum 
relative error which perturbs the coordinates of the points V with respect to V, 
we show, in Appendix B, that 



where Deg(t) is the degree of the term t and du is the maximum degree of the 
terms t\ , . . . , t^- Since both sets V and V are admissible data perturbations, the 
differences between their coordinates is not large and so the value Sr is small. 
We can conclude that, with a suitable choice of the threshold e, the algorithm 
ABM, working on the preprocessed points V, computes, for non-degenerate 
problems, a good approximation to a reference basis. 

Remark. Note that, under small perturbations of V, the matrix Mfc("P) of 
the system @ preserves its full rank, since AT is a reference normal set, but its 
elements slightly change. For this reason, choosing different admissible input 
sets, we obtain sets of polynomials which differ only in their coefficients, in a 
continuous way, and which do not present structural differences. 

5 Extended basis and pseudozeros 

In his book "Numerical Polynomial Algebra" [17] and in several papers (|12j. 
[13] , [14] , [15] , [16] ) Stetter analyzes some aspects of Polynomial Algebra from 
the numerical point of view, when perturbed data and floating point arithmetic 
are used. In this section we compare the set GB with the extended basis devel- 
oped by Stetter. Moreover, we show that, in the terminology of Stetter, V is a 
pseudozero set for the polynomials of GB. 

5.1 The extended basis 

In p~2] , Stetter analyzes the problem of computing the G-basis of an ideal deter- 
mined by a system of polynomials with perturbed coefficients, pointing out that 
the G-basis can be structurally altered by the data uncertainty. The numerical 
instabilities occur, analogously to the case of ideals of points, when the zeros of 
the perturbed polynomials are "near" to a set of points associated to a normal 
set Af p which differs from the perturbed normal set Af p . In order to solve the 
numerical instabilities, Stetter introduces the notion of extended basis GBg, 
that is an approximation to the G-basis, computed using the normal set Af p 



e* < Deg(t)S R and e M < Vk (ImSr, 
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instead of N p . The computed set GBe is not a G-basis in general, but, if there 
exists a G-basis GB structurally corresponding to GBe such that for each poly- 
nomial gE £ GBe there is a polynomial g £ GB with "similar" coefficients, 
then GBe is a "numerically stable" basis of the ideal associated to the input 
polynomials (see theorem 4.1 in [H]). 

Modifying the technique proposed by Stetter, it is possible to compute the 
extended basis for ideals of perturbed points too. Given the set of m points V, 
the normal set Af p coincides with the normal set AT = {Ii,. .. ,t m } computed 
by the ABM algorithm, because AT characterizes all the admissible input sets. 
Following Stetter's method, each polynomial gE £ GBe corresponds to a power 
product t £ C\77] and can be written as 



9E 



t-J^cftj, (6) 



using all the terms of TV". The vector ce of the coefficients of gE is the solution 
of the system 

M m (V)c E = t(V), with M m (V) = [ii(n ■ ■ • ,t m (V)] , (7) 

and so gE vanishes at V . 

Since for each t £ C \N) there exist both a polynomial ~g £ GB and a poly- 
nomial gs £ GBe, there exists an one-to-one correspondence between GB and 
GBe- In addition, since the vector c of the coefficients of ~g is the least square 
solution of the problem ([3]), we have that 



g(V) = t(P) - M k (T)c, and from 0, g(V) = M m (V)c E - M k (P)c. 
Since the matrix Mk(V) consists of the first k columns of M m (7 : '), we have 



g(V) = M m (V)c E - M m (V) 



= M m (P)Ac, where Ac = c E - 



Since Ac = (M^P))- 1 ^) and c E = {M m (T>))~ 1 ~g(P) , we obtain [8] 

«k < || AC || 2 < \mk and \mh < l|c , < p&, 



are respectively the minimum and the maximum singular 
values of Al m ('P). Since the ABM algorithm computes an element of GB if the 
relative residual of the least square problem ([2|) is less than the threshold e, we 
can derive the following upper bound for the differences of the coefficients of 
the elements of GB and of GBe- 

I|AC| ' 2 <MB^^<eK 2 (M m (V)), 



\Ce\\2 \\t{V)\\ 2 O-min 
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where K%{M m (j>)) = o- max /a min is the 2-condition number of M m (V). Since 
in general the previous upper bound overestimates the relative error on the 
coefficients, we can conclude that GB is "near" to the extended basis GBe- 

Remark. The power product t £ C [AT] is probably not the leading term of 
the polynomial gs defined by ©■ In fact, if tk < i < tk+x, since in general 
t(V) is not in the span of the first k columns of M m (V) in the classical sense, 
then the last m — k components of eg are different from zero. In contrast, 
by construction, t is the leading term of g 6 GB and so, with respect to the 
extended basis, GB has the advantage that LT{GB) = C[J7], that is LT(GB) 
coincides with the leading term of a reference basis. 

5.2 Pseudozeros 

Even if sometimes V is a zero set for the polynomials of GB, in general the 
elements of GB do not vanish at V , but we show that the preprocessed points 
are pseudozeros for the polynomials of GB. 

Intuitively, given a tolerance S, a point p is a pseudozero of a polynomial g 
if it is the exact zero of a polynomial g whose coefficients differ from those of g 
by less than S. This idea, formalized by Stetter in [15] . can be generalized for a 
set of points. 

Definition 5.1 Let S be a set of power products. 

A point p is a pseudozero of g = X^ es a 3 xJ w ^ respect to S and tolerance S 
if p is an exact zero of some g £ Ns{g,S), where 

Ns(g,§) = {g ■ g{x) = ajx ,^ -aj\< 5}. 
xies 

A set of point V is a set of pseudozeros of g with respect to S and tolerance 5 if 
it is a set of exact zeros of some g € Ns(g, S). 

A set of point V is a set of pseudozeros for a system of k polynomials <?i, . . . , g^, 
whose supports belong to S , with respect to S and tolerances T> — {Si, . . . ,5k} if 
it is a set of pseudozeros with tolerance 5i for each polynomial g i; i = 1 . . . fc, of 
the system. 

The set V of the preprocessed points is a pseudozero set for GB with respect 
to S = AfUC [AT] , as shown in the following theorem. 

Theorem 5.2 The points of V are pseudozeros of the system consisting 
of the polynomials GB = {g l7 . . . ,g m } with respect to S and with tolerances 
V = {Si, ... , S m }, where 

S=A7UC[A7] and ft = ^i(P)h _ 

^min 
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Proof. Given a polynomial g i = t — Y2j=i c jtj ^ GB, there exists qe — 
t — Yl*jLi c ftj e GBe, vanishing at V, such that 

\cf--c j \<\\MU< M ^ )k = s, 

^rnin 

Since the supports of ~g and §e are contained in S — A/ U C[Af] and gs belongs 
to Ns{gi,Si) we have that V is a pseudozero set of GB. (} 
Each value 5i is small with respect to the norm of the vector c of the coef- 
ficients of <7j. In fact, if g i € GB and gs S G-Bg have similar coefficients, from 
the results of the previous section we have 

Si < eK(M m (V))\\c E \\ and so S t < eK(M m (V))\\c\\ + 0(e 2 ). 

The conclusion follows, since the previous upper bound is, in general, an over- 
estimation and e is a small value of the threshold. 

6 Numerical examples 

In this section we present two numerical examples which illustrate how the 
ABM algorithm approximates a reference basis. From the perturbed points, 
after the preprocessing phase, the ABM algorithm computes in Example 6.1 a 
reference basis and in Example 6.2 an approximation to it. In these examples 
we use the term ordering DegLex; in addition, the coordinates of the points 
and the coefficients of the polynomials are displayed with a finite number of 
digits, but all computations are performed in exact arithmetic using CoCoA 4.2 
[5]. Moreover, in the following examples the threshold in the test of the ABM 
algorithm has been suitably chosen. A more exhaustive survey of numerical 
experiments will be presented in a forthcoming paper. 

Example 6.1 

Given the perturbed points V known with uncertainty sq = 0.1, 

V = {(-2.45, -3.6), (-0.53,-1.45), (1.5,0.45), (3.5,2.5)}, 

the exact DegLex- Grobner basis of the ideal 1(V) is 

( xy - 0.97550y 2 + 1.45450a; - 2.48031y - 1.54307 
GB = < x 2 - 0.95161y 2 + 0.89208a; - 2.94110?/ - 2.07192 

[ y 3 + 0.64549y 2 + 91.31103a; - 98.56564y - 92.83385 

The exact G-basis GB cannot be a reference basis, since there is a small 
relative residual associated to a power product of its normal set Af, as pointed 
out in the following table. 

The small relative residual associated to the power product x means that 
the perturbed points are "almost aligned" and so a normal set corresponding 
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St 


1 y x y 2 


Rel. Res. 


0.97 0.03 0.45 



to the "configuration" of aligned points turns out to be more meaningful and 
stable than N. The ABM algorithm can detect this situation. In fact, using for 
example the threshold e = 0.1 and working on the preprocessed points V 

V = {(-2.475,-3.55), (-0.49,-1.475), (1.475,0.49), (3.55,2.475)}, 

the ABM algorithm computes the set J7 = {1, y, y 2 , y 3 } and the set of polyno- 
mials 

T rs_l x- 0.99979 ? y - 1.02989 

\ y A + 2.06y 3 - 8.45012y 2 - 9.43141y + 6.35026. 

AT" is a reference normal set, since the relative residuals associated to its terms 
are greater than 0.22 and the relative residuals associated to x and y A are less 
than 0.02. In addition, the set GB is a reference basis, since it is the DegLex- 
Grobner basis of the ideal Z("P + ), with normal set 77, where V + is the following 
admissible input set computed by Newton iteration 

V + = {(-2.519, -3.55), (-0.444,-1.475), (1.519, 0.49), (3.504, 2.475)}. 



Example 6.2 

Let V be the set of the perturbed points with uncertainty so = 0.1, 

V= {(0,4.1), (0.05,-4.03), (3.1,0.1), (-3,0.03), 

(2.37,2.5), (-2.4,-2.33), (2.31,-2.486), (-2.4,2.4)}. 

The exact DegLex-Grobner basis of the ideal I(V) is 



GB 



x 2 y 



xy° 



-0.00164a;y 2 + 0.52911y 3 + 0.15329.x 2 - 0.01032x2/ 
+0.05231y 2 - 0.02767a; - 8.75015?/ - 1.47114, 
+0.60405a; 2 / 2 + 0.00120y 3 - 1.38130a; 2 - 0.04581a;?/ 
-0.72231y 2 - 9.17153a; + 0.03536y + 11.91429, 
-2.50589a;y 2 - 0.34459y 3 - 178.88212a; 2 + 1.01537a;y 
-117.12255J/ 2 + 17.87987a; + 11.35992y + 1663.42803, 
-0.09269a;y 2 + 0.04595y 3 + 1.13505a; 2 - 5.90939a;y 
+0.65885y 2 + 0.28526a; - 1.06888y - 9.86020 



The exact normal set N of GB cannot be a reference normal set, since the 
relative residual associated to its power product a; 2 is very small, as reported in 
the following table. 

Since M can change structurally for small data perturbations, GB is not a 
reference basis. The ABM algorithm, using the threshold e = 0.1 and working 
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St 


i y 


2 

x y 


xy 


x y° xy 


Rol. Res. 


0.99 


1 0.64 


0.99 


0.03 0.39 0.57 



on the preprocessed points V 

V= {(0,4.065), (0,-4.065), (3.05,0), (-3.05,0), 

(2.405,2.405), (-2.405,-2.405), (2.405,-2.405), (-2.405,2.405)} 

computes a set of power products 77 — {l,y,x,y 2 ,xy,y 3 ,xy 2 ,y 4 } which is a 
normal set associated to relative residuals greater than 0.22. Since all the rela- 
tive residuals are sufficiently greater or less than the threshold, 77 turns out to 
be a reference normal set. Moreover, the ABM algorithm computes the set of 
polynomials 

f x 2 + 0.55840y 2 - 9.13935, 

GB = < xy 3 - 5.78402x7/, 

[ y 5 - 22.30825y 3 + 95.57653y, 

whose zero set {(3.02313, 0), (—3.02313, 0)} is not an admissible input set. Nev- 
ertheless, Newton method applied to the function F, defined by follows 
descent directions for For this reason the Newton iteration computes the 

admissible input set 

V + = {(0,4.065), (0,-4.065), (3.02,0), (-3.02,0), 

(2.41,2.41), (-2.41,-2.41), (2.41,-2.41), (-2.41,2.41)}. 

This fact suggests that the threshold has been suitably chosen and that 77 
and GB are a good approximation to the exact case. In fact, the set GB 
approximates the set GB 

( x 2 + 0.5625y 2 - 9, 
GB = I xy 3 - 5.76x2/, 

[ y 5 - 21.76y 3 + 92.16y 

which is a reference basis, since it is a DegLex-Grobner basis, with normal set 
77, of the ideal of the admissible input points 

{(0, 4), (0, -4), (3, 0), (-3, 0), (2.4, 2.4), (-2.4, -2.4), (2.4, -2.4), (-2.4, 2.4)}. 



Appendix A: the preprocessing phase 

In this appendix we summarize the preprocessing phase of the perturbed points 
V . This procedure finds the coordinates which differ by less than the error 
estimation and replaces them by a single representative value. For this reason 
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the set V of the preprocessed points is equivalent to V from the computational 
point of view, even if it can happen that the cardinality of V is less than the 
cardinality of V. In order to preserve or to get back the "geometrical symme- 
tries" of the exact points, the preprocessing algorithm works on the absolute 
values of the coordinates of all the points, assembled together in a set Y, and it 
computes a preprocessed set Y. The preprocessed points V are built from the 
set Y in an obvious way. 

Let So be an estimate of the errors on the data and let Y = {y\, . . . ,y n } be 
the set of the absolute values of the coordinates of the input points, sorted in 
non-decreasing order. If (k + 1) elements yj, . . . ,yj+k of Y are such that the 
intersection 

j+k 

is not empty, then any element y~j £ Xj can represent, from the computational 
point of view, the values yj, . . . , yj+k, since it differs from them by less than sq. 
For this reason we replace yj, . . . , yj+k with a value ljj £ Xj in the set Y of the 
preprocessed coordinates. 

The preprocessing algorithm for building the set Y can be described as 
follows. First of all, each element of Y less than so is removed from Y and 
the value is inserted into Y. After this step, the algorithm computes for each 
element yj £ Y the largest non empty intersection Xj, and then it processes the 
set Xj which contains the maximum number of elements of Y, that is the set 
with the maximum index k. If there are several intersections with this property, 
the set with the minimum index j is chosen. If Xj is the intersection to process, 
all the elements yj, . . . , yj+k <= Xj are removed from Y and the middle point y^ 
of Xj is inserted in Y. 

The preprocessing algorithm ends when Y is empty or when it contains 
elements which differ by more than so. In this case such values do not require 
any preprocessing treatment and so they are removed from Y and inserted 
directly into Y. 

Even though it can happen that [y, — s , y l + s ] n [y i+1 — s , y i+1 + s ] 7^ 0, 
we decided not to repeat the preprocessing phase on the set Y, since all the new 
coordinates differ by more than so- In fact, if the set Xj° is processed, 

[yj-! - s , yj-! + s ] n X* = 0, [yj+k+i - s , y 3 +k+i + s ] n Xf = 0, 

and so the middle point y^ of Xj differs from yj-! and yj+k+i by more than 
so. Since yj+k+i — Vj-i > 2so, the preprocessing phase works separately on 
the sets Yi = {j/i, . . . , j/j-i} and Y2 = {yj+k+! , • • • , y n }- The values obtained 
by preprocessing Yi and F 2 are, respectively, less than the maximum of Y\ and 
greater than the minimum of Y 2 , and so they differ from y~j by more than s . 
Analogously, analyzing the sets Yi and Yi separately, we can show that all the 
preprocessed values differ to each other by more than sq, that is they are "well 
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separated" . This result points out that it is computationally better to use the 
set V than the set V as input points. 

Appendix B: an upper bound for e t and cm 

In this section we present an upper bound for the values e t and cm, useful for 
estimating how well the set GB approximates a reference basis. 

First of all, we analyze the sensitivity of a term evaluated at a perturbed 
point. Let t = or} 1 . . . x{° be a power product belonging to R[xi . . . x s ] and let 
Deg(t) = ji + . . . + is be its degree. If p = (p x , . . . ,p s ) is a perturbation of the 
point p — (pi, . . . ,p s ) such that p i = p,(l + Si), |<5,| < 5 P , i = 1 . . . s, then 

t(p) = pf . . . j%> = p] 1 (1 + fc)* . . . pj- (1 + 6 a y- = t(p)(l + 5^ . . . (1 + 5.Y' . 

By a first-order error analysis, ignoring the errors of higher order, we obtain 

\t(p) - t( P )\ < |t(p)|(ji|Ji| + • ..js\Ss\) < \t(p)\Deg(t)5 p . 

Now, we can easily upper bound e t as follows. Let V = {p 1 , . . . ,p m } be a 
set of points perturbation of V = {pi, . . . ,p m }, such that pi — (jPi\ ■ ■ ■ ,pi^) 
and p t = (pW(i + s[ l) ), . . . ,pW(l + #)), where |5j*>| < 5 R , i = I . . .m and 
j = 1. ..s. Then we have 

m 

\\t(V)-t(V)f 2 <Yj( Pl fDeg(tf5l < Deg(t) 2 S 2 R \\t(V)\\ 2 that is e t < Deg(t)S R . 

Analogously, we can also upper bound cm- Let M(T) and M(P) be the matrices 
whose columns are generated by the terms t\ . . . tk evaluated, respectively, at 
V and V. If du = ma,x{Deg(ti), . . . , Deg(tk)} is the maximum degree of the 
terms t\ . . . tk, we obtain 

\\M(ryM(r)\\% = £ ll^-^lli < 4^ 2 m £ H*^)Hi - <P m 8 r \\m{v)\\ 2 f , 

3=1 3=1 

where || \\f is the Frobenius matrix norm. Since, given an m x fc matrix A, we 
have ||A|| 2 < \\A\\ F < Vk\\A\\ 2 , then 

||M(?) - M(V)\\ 2 , /r WMCP) - M(V)\\ F rr , 
£M = \\M(V) h ~ ^ \\M(V)\\ F * VkdM6R - 
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